home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48_1 / eign.stu < prev    next >
Text File  |  1995-03-23  |  7KB  |  409 lines
  1. Article 4537 of comp.sys.handhelds:
  2. Newsgroups: comp.sys.handhelds
  3. Path: en.ecn.purdue.edu!dominop
  4. From: dominop@en.ecn.purdue.edu (Philippos A. Peleties)
  5. Subject: Re: Eigen-values/vectors request
  6. Message-ID: <1991Feb25.145302.9837@en.ecn.purdue.edu>
  7. Keywords: Eigenvalues, Schur form
  8. Organization: Purdue University Engineering Computer Network
  9. References: <1991Feb24.115734.14198@cat.de>
  10. Date: Mon, 25 Feb 91 14:53:02 GMT
  11.  
  12. Well, here are my versions of eigenvalue routines.
  13.  
  14. The first one is through the Leverierre algorithm. Not very good
  15. (roots of polynomials are sensitive to coefficient perturbations) but 
  16. quite adequate for quick checks.
  17.  
  18. The second routine converts a matrix to a schur form. Use 0 or 1 to
  19. indicate whether you want shift or not. Shift results in quadratic
  20. convergence rate. Sometimes you might not want to use shift though.
  21.  
  22. As far as eigenvectors are concerned, that's not so easy. I don't
  23. have any eigenvector routines ( null(A-lambda*I) is not the way to
  24. go). Does anybody have any intelligent algorithms he/she wants to
  25. provide?
  26.  
  27.  
  28.  
  29. Have fun.
  30.  
  31. ----------------------------------------------------------------------
  32.  
  33. Name:        EIGEN
  34. ----
  35.  
  36. Status:        In hp28s.
  37. ------
  38.  
  39. Author:        Philip Peleties.
  40. ------
  41.  
  42. Function:    Compute the eigenvalues of a matrix via the
  43. --------    Leverierre algorithm.
  44.  
  45. Input:        A    Matrix.
  46. -----
  47.  
  48. Output:        List with eigenvalues.
  49. ------
  50.  
  51. Calls:        TRACE, BAIR
  52. -----
  53.  
  54. Checksum:    DFBF
  55. --------
  56.  
  57. << -> A
  58.    << A SIZE LIST->
  59. DROP IDN -> n N
  60.       << 1 1 n
  61.          FOR i N A *
  62. TRACE NEG i / DUP n
  63. IDN * N A * + 'N'
  64. STO
  65.       NEXT 1 n +
  66. ->LIST BAIR
  67.       >>
  68.    >>
  69. >>
  70.  
  71.  
  72. ------------------------------------------------------------------
  73.  
  74.  
  75. Name:        TRACE
  76. ----
  77.  
  78. Status:        In hp28s.
  79. ------
  80.  
  81. Author:        Philip Peleties.
  82. ------
  83.  
  84. Function:    Computes the trace of a square matrix.
  85. --------
  86.  
  87. Input:        A    Square matrix.
  88. -----
  89.  
  90. Output:        Trace of A.
  91. ------
  92.  
  93. Calls:        None.
  94. -----
  95.  
  96. Checksum:    A803
  97. --------
  98.  
  99. << -> A
  100.    << 0 A SIZE LIST->
  101. DROP2 1 SWAP
  102.     FOR i A i i 2
  103. ->LIST GET +
  104.     NEXT
  105. >>
  106.  
  107.  
  108. --------------------------------------------------------------------
  109.  
  110. Name:        BAIR
  111. ----
  112.  
  113. Status:        In hp28s.
  114. ------
  115.  
  116. Author:        Mark Wicks.
  117. ------
  118.  
  119. Function:    Computes the roots, both real and complex, of a 
  120. --------    polynomial.
  121.  
  122. Input:        List containing coefficients of polynomial.
  123. -----
  124.  
  125. Output:        List containing roots of the polynomial.
  126. ------
  127.  
  128. Calls:        NR.
  129. -----
  130.  
  131. Checksum:    9AF7.
  132. --------
  133.  
  134. << { } SWAP 0 0 0 0 0
  135. 0 0 0 0 0 0
  136. [ 0 0 ]
  137. -> a b r s p q b1 b2
  138. bp1 bp2 bq1 bq2 opq
  139.   <<
  140.     IF a SIZE DUP 2
  141. / IP 2 * ==
  142.     THEN a NR SWAP
  143. 'a' STO +
  144.     END
  145.     WHILE a SIZE 3 >
  146.     REPEAT 'a' 2 GET
  147. 'a' 1 GET / 'p' STO
  148. 'a' 3 GET 'a' 1 GET
  149. / 'q' STO 
  150.       IF p 0 ==
  151.       THEN .1 'p'
  152. STO
  153.       END
  154.       IF q 0 ==
  155.       THEN .2 'q'
  156. STO
  157.       END
  158.       DO 'a' 1 GET
  159. 'b2' STO 'a' 2 GET p
  160. b2 * - 'b1' STO 0
  161. 'bp2' STO 0 'bq2'
  162. STO b2 NEG 'bp1' STO
  163. 0 'bq1' STO b2 b1 2
  164. ->LIST 'b' STO 3 a
  165. SIZE 2 - DUP2
  166.     IF <=
  167.         THEN
  168.       FOR i 'a'
  169. i GET p b1 * - q b2
  170. * - b1 NEG p bp1 * -
  171. q bp2 * - bp1 'bp2'
  172. STO 'bp1' STO b2 NEG
  173. p bq1 * - q bq2 * -
  174. bq1 'bq2' STO 'bq1'
  175. STO b1 'b2' STO 'b1'
  176. STO b b1 + 'b' STO
  177.       NEXT
  178.         ELSE DROP2
  179.     END a DUP
  180. SIZE 1 - GET p b1 *
  181. - q b2 * - 'r' STO a
  182. DUP SIZE GET q b1 *
  183. - 's' STO p q 2
  184. ->ARRY r s 2 ->ARRY b1
  185. NEG p bp1 * - q bp2
  186. * - b2 NEG p bq1 * -
  187. q bq2 * - q NEG bp1
  188. * b1 NEG q bq1 * - {
  189. 2 2 } ->ARRY / -
  190. ARRY-> DROP 'q' STO
  191. 'p' STO
  192.       UNTIL p q 2
  193. ->ARRY DUP opq - ABS
  194. SWAP DUP 'opq' STO
  195. ABS .0000000001 * <
  196.       END p NEG 2 /
  197. DUP SQ q - $sqrt$ DUP2 +
  198. 3 ROLLD - 2 ->LIST +
  199. b 'a' STO
  200.     END 'a' 2 GET
  201. NEG 'a' 1 GET / 2 / 
  202. DUP SQ 'a' 3 GET 'a'
  203. 1 GET / - $sqrt$ DUP2 + 3
  204. ROLLD - 2 ->LIST +
  205.   >>
  206. >>
  207.  
  208.  
  209. --------------------------------------------------------------------
  210.  
  211. Name:        NR
  212. ----
  213.  
  214. Status:        In hp28s.
  215. ------
  216.  
  217. Author:        Mark Wicks.
  218. ------
  219.  
  220. Function:    Computes a single real root of a polynomial and returns 
  221. --------    the deflated polynomial and the root.
  222.  
  223. Input:        List containing coefficients of polynomial.
  224. -----
  225.  
  226. Output:        A list containing the coefficients of the deflated 
  227. ------        polynomial.
  228.         A single real root of the polynomial.
  229.  
  230. Calls:        None.
  231. -----
  232.  
  233. Checksum:    B7C2
  234. --------
  235.  
  236. << 0 0 0 0 0 -> a c r
  237.  p b db
  238.   << 'a' 2 GET 'a' 1
  239. GET / 'p' STO
  240.     DO 0 'db' STO a
  241. 1 GET 'b' STO b 1
  242. ->LIST 2 a SIZE 1 -
  243.       FOR i 'a' i
  244. GET p b * - db p * b
  245. + NEG 'db' STO DUP
  246. 'b' STO +
  247.       NEXT 'c' STO
  248. 'a' a SIZE GET p b *
  249. - b p db * + NEG / p
  250. SWAP DUP2 - 'p' STO
  251.     UNTIL ABS SWAP
  252. ABS .0000000001 * <
  253.     END c p NEG
  254.   >>
  255. >>
  256.  
  257. ----------------------------------------------------------------------
  258.  
  259. Name:        SCHUR
  260. ----
  261.  
  262. Status:        In hp28s.
  263. ------
  264.  
  265. Author:        Philip Peleties.
  266. ------
  267.  
  268. Function:    Compute the Schur form of a matrix.
  269. --------
  270.  
  271. Input:        A    Matrix.
  272. -----
  273.         r    Number of QR iterations.
  274.         s    Shift indicator (0 if no shift, 1 if shift).
  275.  
  276. Output:        Q    Orthogonal matrix.
  277. ------
  278.         T    Schur form (A=Q*T*Q')
  279.                 number    ABS(Q*T*Q'-A)
  280.  
  281. Calls:        QR
  282. -----
  283.  
  284. Checksum:    5E5E
  285. --------
  286.  
  287.  
  288. << -> A r s
  289.    << A SIZE LIST->
  290. DROP IDN DUP
  291.       IF s 0 ==
  292.       THEN 0 *
  293.       END A -> n Q IDNT
  294. T
  295.    << 1 r
  296.       START Q T 'T(n
  297. ,n)' EVAL IDNT * -
  298. QR DROP DUP DUP TRN
  299. T ROT * * 'T' STO *
  300. 'Q' STO
  301.       NEXT Q T Q T Q 
  302. TRN * * A - ABS
  303.     >>
  304.   >>
  305. >>
  306.  
  307. ----------------------------------------------------------------
  308.  
  309. Name:        QR
  310. ----
  311.  
  312. Status:        In hp28s.
  313. ------
  314.  
  315. Author:        Philip Peleties.
  316. ------
  317.  
  318. Function:    Compute the QR factorization of a matrix.
  319. --------
  320.  
  321. Input:        A    Matrix.
  322. -----        
  323.  
  324. Output:        Q    Orthogonal matrix.
  325. ------        R    Upper triangular matrix (A=Q*R).
  326.  
  327. Calls:        None.
  328. -----
  329.  
  330. Checksum:    7937
  331. --------
  332.  
  333. << -> A
  334.    << A SIZE LIST->
  335. DROP -> n m
  336.      << n IDN DUP A 
  337. TYPE
  338.         IF 4 ==
  339.         THEN (1,0) *
  340.         END A 0 0 -> Q 
  341. IDNT R s c
  342.        << 1 m
  343.          FOR j 1 j +
  344. n DUP2
  345.            IF >
  346.            THEN DROP2
  347.            ELSE
  348.              FOR i 'R
  349. (j,j)' EVAL 'R(i,j)' 
  350. EVAL DUP ABS
  351.               IF 0 
  352. ==
  353.               THEN 1
  354. 'c' STO 0 's' STO
  355. DROP2
  356.               ELSE
  357. DUP2 ABS SWAP ABS 
  358. IF >
  359. THEN ABS / DUP ABS
  360. SQ 1 + $sqrt$ INV DUP 'R(
  361. i,j)' EVAL DUP ABS /
  362. * 's' STO * 'c' STO
  363. ELSE SWAP ABS / DUP
  364. ABS SQ 1 + $sqrt$ INV DUP
  365. 'R(j,j)' EVAL DUP
  366. ABS / * 'c' STO * 
  367. 's' STO
  368. END
  369.               END Q
  370. IDNT i i 2 ->LIST c 
  371. PUT i j 2 ->LIST s
  372. NEG PUT j j 2 ->LIST 
  373. c CONJ PUT j i 2 
  374. ->LIST s CONJ PUT TRN 
  375. DUP TRN R * 'R' STO 
  376. * 'Q' STO
  377.             NEXT
  378.           END
  379.         NEXT Q R
  380.         IF TYPE 4 ==
  381.         THEN IDNT n 
  382. m MIN DUP 2 ->LIST 
  383. DUP R SWAP GET DUP 
  384. ABS
  385.          IF
  386. .000000000001 <=
  387.          THEN 3
  388. DROPN
  389.          ELSE DUP
  390. ABS / PUT DUP 3 
  391. ROLLD * SWAP TRN R *
  392. 'R' STO
  393.          END
  394.        END R
  395.       >>
  396.     >>
  397.   >>
  398. >>
  399.  
  400.  
  401. ---------------------------------------------------------------------------------
  402. -- 
  403.  
  404.  
  405. I speak for myself, I think for myself, I work for myself ... but I don't want
  406. to play by myself ... so bring your toys and let's share ...
  407.  
  408.  
  409.